dsMERGIATO <- read.csv("/Users/leomarcellopoli/Documents/Nonparametric/Project/Non-parametric-Statistics-2025-2026/dataset/dataset_MERGIATO", header = TRUE, sep = ",")
{#codice per trasformare data in numerico
dsMERGIATO <- dsMERGIATO[dsMERGIATO$Lat >= 45, ]
dsMERGIATO <- dsMERGIATO %>%
  mutate(
    Date = as.Date(Date),
    Date_numeric = as.numeric(Date - as.Date("2021-12-31")) # as.Date("2021-12-31") è lo 0
  ) %>%
  select(Date_numeric, everything(), -Date)
}
time_grid <- sort(unique(dsMERGIATO$Date_numeric))
coordinates <- unique(group_by(dsMERGIATO[,2:3]))
vars_to_func <- c("Chl", "Temp", "Salinity", "MLD", "uo", "vo", "Current_Speed", "Solar_Flux", "Heat_Flux")
coord_reference <- dsMERGIATO %>%
  distinct(Lon, Lat) %>%
  arrange(Lon, Lat) # Deve essere identico all'arrange usato nel pivot_wider!
world_map <- map_data("world", region = c("Italy", "Croatia", "Slovenia"))
temporary_ds <- dsMERGIATO %>%
  filter(Lon == coordinates$Lon[1] & Lat == coordinates$Lat[1]) %>%
  arrange(Date_numeric) %>%
  select(Date_numeric, all_of(vars_to_func[1]))
fdata_Chl_FIRST_COORDINATE <- fData(grid = time_grid, values = matrix(temporary_ds$Chl, nrow = 1))
plot(fdata_Chl_FIRST_COORDINATE, main = "Chl at first coordinate")

Removing the useless variables:

Ora lo facciamo per tutte le coordinate e tutte le variabili (con lo facciamo, intendo dire io e Gemini).

Plot of variables during the year

for(var_name in vars_to_func) {
  plot(dsfunctional[[var_name]], main = paste("Tutte le curve funzionali di", var_name), ylab = var_name)
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
}

Divisione zone mare

Plot like this mean almost nothing, it could make sense to do clustering and to indentify the coordinates relative to bathing areas. Anyway I will to some clustering in a deterministic way, by selecting the coords based on the values of the chlorophyll.

Depth measures

Chlorophyll

For now we’ll work only on chlorophyll, but we could extend the analysis to other variables.

fdata_Chl <- dsfunctional[["Chl"]]
band_depth_Chl <- BD(fdata_Chl)
modified_band_depth_Chl <- MBD(fdata_Chl)
which.max(band_depth_Chl)
## [1] 36
which.min(band_depth_Chl)
## [1] 1
idx_max_bd <- which.max(band_depth_Chl)
idx_min_bd <- which.min(band_depth_Chl)
coord_max_bd <- coord_reference[idx_max_bd, ]
coord_min_bd <- coord_reference[idx_min_bd, ]
print(paste("Max BD (Curva più centrale/tipica) si trova a:", coord_max_bd$Lat, coord_max_bd$Lon))
## [1] "Max BD (Curva più centrale/tipica) si trova a: 45.229 12.5"
print(paste("Min BD (Curva più atipica/outlier) si trova a:", coord_min_bd$Lat, coord_min_bd$Lon))
## [1] "Min BD (Curva più atipica/outlier) si trova a: 45.188 12.333"
{
  # Plot di base con tutte le curve
  plot(fdata_Chl, 
       main = "Curva funzionale di Chl: Max (Blu) e Min (Rosso)", 
       ylab = "Chl", 
       ylim = c(0, max(fdata_Chl$values)),
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_Chl[idx_max_bd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_Chl[idx_min_bd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max BD (Tipico)", "Min BD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

which.max(modified_band_depth_Chl)
## [1] 136
which.min(modified_band_depth_Chl)
## [1] 1
idx_max_mbd <- which.max(modified_band_depth_Chl)
idx_min_mbd <- which.min(modified_band_depth_Chl)
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]
print(paste("Max MBD si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD si trova a: 45.396 12.833"
print(paste("Min MBD si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD si trova a: 45.188 12.333"
{
  # Plot di base con tutte le curve
  plot(fdata_Chl, 
       main = "Curva funzionale di Chl: Max (Blu) e Min (Rosso)", 
       ylab = "Chl", 
       ylim = c(0, max(fdata_Chl$values)),
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_Chl[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_Chl[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Let’s see this coordinates on the map:

# Plot veloce
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max Depth (Il punto più "normale") in BLU
points(coord_max_bd$Lon, coord_max_bd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_bd$Lon, coord_max_bd$Lat, labels = "Max BD", pos = 3, col="blue")

# Aggiungi Min Depth (L'Outlier) in ROSSO
points(coord_min_bd$Lon, coord_min_bd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_bd$Lon, coord_min_bd$Lat, labels = "Min BD", pos = 3, col="red")

title("Posizione della Curva Mediana (Blu) e dell'Outlier (Rosso)")

# Plot veloce
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max Depth (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min Depth (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("Posizione della Curva Mediana (Blu) e dell'Outlier (Rosso)")

Map of all band depth

data_plot_bd <- coord_reference %>%
  mutate(Depth_Value = band_depth_Chl)

ggplot() +
  # Sfondo grigio (Italia)
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla Depth
  geom_point(data = data_plot_bd, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori: Viridis è ottima (Giallo=Alto/Centrale, Viola=Basso/Outlier)
  scale_color_viridis_c(option = "viridis", name = "Band Depth") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Band Depth (BD) - Clorofilla",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Map of all modified band depth of Chloophyll

data_plot_mbd <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_Chl)

ggplot() +
  # Sfondo grigio (Italia)
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla Depth
  geom_point(data = data_plot_mbd, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori: Viridis è ottima (Giallo=Alto/Centrale, Viola=Basso/Outlier)
  scale_color_viridis_c(option = "viridis", name = "Modified Band Depth") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - Clorofilla",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Temperature

Now we’ll switch to analysing Temperature

fdata_Temp <- dsfunctional[["Temp"]]
modified_band_depth_Temp <- MBD(fdata_Temp)
# Identificazione indici Max e Min
idx_max_mbd <- which.max(modified_band_depth_Temp)
idx_min_mbd <- which.min(modified_band_depth_Temp)

# Coordinate corrispondenti
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]

# Print coordinate
print(paste("Max MBD (Temp) si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD (Temp) si trova a: 45.396 13.042"
print(paste("Min MBD (Temp) si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD (Temp) si trova a: 45.396 12.375"
# Plot delle curve
# Uso range() per ylim per adattarmi alla scala della temperatura
y_limits <- range(fdata_Temp$values) 

{
  # Plot di base con tutte le curve
  plot(fdata_Temp, 
       main = "Curva funzionale di Temp: Max (Blu) e Min (Rosso)", 
       ylab = "Temp", 
       ylim = c(0, max(fdata_Temp$values)),
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_Temp[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_Temp[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Map of all modified band depth for Temperature

# Plot veloce statico
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max MBD (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min MBD (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("Temp: Posizione Curva Mediana (Blu) e Outlier (Rosso)")

Map of all modified band depth for Temperature

# Creazione dataframe per ggplot specifico per Temp
data_plot_mbd_temp <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_Temp)

ggplot() +
  # Sfondo grigio (Italia) - usa world_map già caricato in precedenza
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla MBD della Temperatura
  geom_point(data = data_plot_mbd_temp, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori Viridis
  scale_color_viridis_c(option = "viridis", name = "MBD Temp") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - Temperatura",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Salinity

Now we’ll switch to analysing Salinity

fdata_Salinity <- dsfunctional[["Salinity"]]
modified_band_depth_Salinity <- MBD(fdata_Salinity)
# Identificazione indici Max e Min
idx_max_mbd <- which.max(modified_band_depth_Salinity)
idx_min_mbd <- which.min(modified_band_depth_Salinity)

# Coordinate corrispondenti
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]

# Print coordinate
print(paste("Max MBD (Salinity) si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD (Salinity) si trova a: 45.563 13.208"
print(paste("Min MBD (Salinity) si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD (Salinity) si trova a: 45.188 12.333"
# Plot delle curve
# Uso range() per ylim per adattarmi alla scala della salinità
y_limits <- range(fdata_Salinity$values) 

{
  # Plot di base con tutte le curve
  plot(fdata_Salinity, 
       main = "Curva funzionale di Salinity: Max (Blu) e Min (Rosso)", 
       ylab = "Salinity", 
       ylim = c(min(fdata_Salinity$values), max(fdata_Salinity$values)),
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_Salinity[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_Salinity[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Map of all modified band depth for Salinity

# Plot veloce statico
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max MBD (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min MBD (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("Salinity: Posizione Curva Mediana (Blu) e Outlier (Rosso)")

Map of all modified band depth for Salinity

# Creazione dataframe per ggplot specifico per Salinity
data_plot_mbd_salinity <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_Salinity)

ggplot() +
  # Sfondo grigio (Italia) - usa world_map già caricato in precedenza
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla MBD della Salinità
  geom_point(data = data_plot_mbd_salinity, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori Viridis
  scale_color_viridis_c(option = "viridis", name = "MBD Salinity") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - Salinity",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Mixed Layer Depth

Now we’ll switch to analysing MLD

# Estrazione dati funzionali MLD
fdata_MLD <- dsfunctional[["MLD"]]

# Calcolo SOLO Modified Band Depth (MBD)
modified_band_depth_MLD <- MBD(fdata_MLD)
# Identificazione indici Max e Min
idx_max_mbd <- which.max(modified_band_depth_MLD)
idx_min_mbd <- which.min(modified_band_depth_MLD)

# Coordinate corrispondenti
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]

# Print coordinate
print(paste("Max MBD (MLD) si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD (MLD) si trova a: 45.563 13.5"
print(paste("Min MBD (MLD) si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD (MLD) si trova a: 45.646 13.125"
# Plot delle curve
# Uso range() per ylim per adattarmi alla scala della MLD
y_limits <- range(fdata_MLD$values) 

{
  # Plot di base con tutte le curve
  plot(fdata_MLD, 
       main = "Curva funzionale di MLD: Max (Blu) e Min (Rosso)", 
       ylab = "MLD", 
       ylim = c(0, max(fdata_MLD$values)),
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_MLD[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_MLD[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Map of all modified band depth for MLD

# Plot veloce statico
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max MBD (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min MBD (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("MLD: Posizione Curva Mediana (Blu) e Outlier (Rosso)")

Map of all modified band depth for MLD

# Creazione dataframe per ggplot specifico per MLD
data_plot_mbd_mld <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_MLD)

ggplot() +
  # Sfondo grigio (Italia) - usa world_map già caricato in precedenza
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla MBD della MLD
  geom_point(data = data_plot_mbd_mld, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori Viridis
  scale_color_viridis_c(option = "viridis", name = "MBD MLD") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - MLD",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Current Speed

Now we’ll switch to analysing Current Speed

# Estrazione dati funzionali Current_Speed
fdata_CS <- dsfunctional[["Current_Speed"]]

# Calcolo SOLO Modified Band Depth (MBD)
modified_band_depth_CS <- MBD(fdata_CS)
# Identificazione indici Max e Min
idx_max_mbd <- which.max(modified_band_depth_CS)
idx_min_mbd <- which.min(modified_band_depth_CS)

# Coordinate corrispondenti
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]

# Print coordinate
print(paste("Max MBD (Current Speed) si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD (Current Speed) si trova a: 45.229 13.292"
print(paste("Min MBD (Current Speed) si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD (Current Speed) si trova a: 45.021 12.625"
# Plot delle curve
# Uso range() per essere sicuro di coprire min e max corretti
y_limits <- range(fdata_CS$values) 

{
  # Plot di base con tutte le curve
  plot(fdata_CS, 
       main = "Curva funzionale di Current Speed: Max (Blu) e Min (Rosso)", 
       ylab = "Current Speed", 
       ylim = y_limits,
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_CS[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_CS[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Map of all modified band depth for Current Speed

# Plot veloce statico
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max MBD (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min MBD (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("Current Speed: Posizione Curva Mediana (Blu) e Outlier (Rosso)")

Map of all modified band depth for Current Speed

# Creazione dataframe per ggplot specifico per Current Speed
data_plot_mbd_cs <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_CS)

ggplot() +
  # Sfondo grigio (Italia) - usa world_map già caricato in precedenza
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla MBD della Current Speed
  geom_point(data = data_plot_mbd_cs, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori Viridis
  scale_color_viridis_c(option = "viridis", name = "MBD Current Speed") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - Current Speed",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Solar Flux

Now we’ll switch to analysing Solar Flux

# Estrazione dati funzionali Solar_Flux
fdata_SF <- dsfunctional[["Solar_Flux"]]

# Calcolo SOLO Modified Band Depth (MBD)
modified_band_depth_SF <- MBD(fdata_SF)
# Identificazione indici Max e Min
idx_max_mbd <- which.max(modified_band_depth_SF)
idx_min_mbd <- which.min(modified_band_depth_SF)

# Coordinate corrispondenti
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]

# Print coordinate
print(paste("Max MBD (Solar Flux) si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD (Solar Flux) si trova a: 45.313 13.042"
print(paste("Min MBD (Solar Flux) si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD (Solar Flux) si trova a: 45.729 13.542"
# Plot delle curve
# Uso range() per essere sicuro di coprire min e max corretti
y_limits <- range(fdata_SF$values) 

{
  # Plot di base con tutte le curve
  plot(fdata_SF, 
       main = "Curva funzionale di Solar Flux: Max (Blu) e Min (Rosso)", 
       ylab = "Solar Flux", 
       ylim = y_limits,
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_SF[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_SF[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Map of all modified band depth for Solar Flux

# Plot veloce statico
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max MBD (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min MBD (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("Solar Flux: Posizione Curva Mediana (Blu) e Outlier (Rosso)")

Map of all modified band depth for Solar Flux

# Creazione dataframe per ggplot specifico per Solar Flux
data_plot_mbd_sf <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_SF)

ggplot() +
  # Sfondo grigio (Italia) - usa world_map già caricato in precedenza
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla MBD del Solar Flux
  geom_point(data = data_plot_mbd_sf, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori Viridis
  scale_color_viridis_c(option = "viridis", name = "MBD Solar Flux") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - Solar Flux",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Heat Flux

Now we’ll switch to analysing Heat Flux

# Estrazione dati funzionali Heat_Flux
fdata_HF <- dsfunctional[["Heat_Flux"]]

# Calcolo SOLO Modified Band Depth (MBD)
modified_band_depth_HF <- MBD(fdata_HF)
# Identificazione indici Max e Min
idx_max_mbd <- which.max(modified_band_depth_HF)
idx_min_mbd <- which.min(modified_band_depth_HF)

# Coordinate corrispondenti
coord_max_mbd <- coord_reference[idx_max_mbd, ]
coord_min_mbd <- coord_reference[idx_min_mbd, ]

# Print coordinate
print(paste("Max MBD (Heat Flux) si trova a:", coord_max_mbd$Lat, coord_max_mbd$Lon))
## [1] "Max MBD (Heat Flux) si trova a: 45.396 13.083"
print(paste("Min MBD (Heat Flux) si trova a:", coord_min_mbd$Lat, coord_min_mbd$Lon))
## [1] "Min MBD (Heat Flux) si trova a: 45.146 12.375"
# Plot delle curve
# Uso range() per essere sicuro di coprire min e max (Heat Flux può essere negativo)
y_limits <- range(fdata_HF$values) 

{
  # Plot di base con tutte le curve
  plot(fdata_HF, 
       main = "Curva funzionale di Heat Flux: Max (Blu) e Min (Rosso)", 
       ylab = "Heat Flux", 
       ylim = y_limits,
       col = "gray80") 
  
  # linea del MAX Band Depth (BLU)
  plot(fdata_HF[idx_max_mbd], col = "blue", lwd = 3, add = TRUE)  
  
  # linea del MIN Band Depth (ROSSO)
  plot(fdata_HF[idx_min_mbd], col = "red", lwd = 3, add = TRUE)
  
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 9),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata)
         lwd = 0.5)                  # Spessore sottile
  
  legend("topright", legend = c("Max MBD (Tipico)", "Min MBD (Outlier)"),
         col = c("blue", "red"), lwd = 3, bg = "white")
}

Map of all modified band depth for Heat Flux

# Plot veloce statico
map('world', xlim = c(12, 14.5), ylim = c(45, 46), col = "gray")
map.axes()

# Aggiungi Max MBD (Il punto più "normale") in BLU
points(coord_max_mbd$Lon, coord_max_mbd$Lat, col = "blue", pch = 19, cex = 2)
text(coord_max_mbd$Lon, coord_max_mbd$Lat, labels = "Max MBD", pos = 3, col="blue")

# Aggiungi Min MBD (L'Outlier) in ROSSO
points(coord_min_mbd$Lon, coord_min_mbd$Lat, col = "red", pch = 19, cex = 2)
text(coord_min_mbd$Lon, coord_min_mbd$Lat, labels = "Min MBD", pos = 3, col="red")

title("Heat Flux: Posizione Curva Mediana (Blu) e Outlier (Rosso)")

Map of all modified band depth for Heat Flux

# Creazione dataframe per ggplot specifico per Heat Flux
data_plot_mbd_hf <- coord_reference %>%
  mutate(Depth_Value = modified_band_depth_HF)

ggplot() +
  # Sfondo grigio (Italia) - usa world_map già caricato in precedenza
  geom_polygon(data = world_map, aes(x = long, y = lat, group = group), 
               fill = "lightgray", color = "white") +
  
  # Punti colorati in base alla MBD del Heat Flux
  geom_point(data = data_plot_mbd_hf, aes(x = Lon, y = Lat, color = Depth_Value), 
             size = 2, alpha = 0.8) +
  
  # Scala colori Viridis
  scale_color_viridis_c(option = "viridis", name = "MBD Heat Flux") +
  
  # Focus sull'Alto Adriatico
  coord_fixed(xlim = c(12, 14.5), ylim = c(45, 46), ratio = 1.3) +
  
  labs(title = "Mappa Modified Band Depth (MBD) - Heat Flux",
       subtitle = "Giallo = Comportamento Tipico | Viola = Comportamento Anomalo",
       x = "Longitudine", y = "Latitudine") +
  theme_minimal()

Functional Boxplot

Finally, we can create functional boxplots for each variable to visualize the distribution of functional data.

for(var_name in vars_to_func) {
  invisible(fbplot(dsfunctional[[var_name]], main = paste("Boxplot di", var_name), ylab = var_name))
  # 2. Aggiunta linee verticali (Griglia temporale)
  abline(v = seq(151, 516, by = 7),  # Sequenza da 151 a 516 ogni 9
         col = "gray",               # Colore grigio
         lty = 2,                    # Tipo di linea (2 = tratteggiata/dashed)
         lwd = 0.5)                  # Spessore sottile
}